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Abstract 

In this project^, we first study the Gaussian-based hidden 
Markov random field (HMRF) model and its expectation- 
maximization (EM) algorithm. Then we generalize it to 
Gaussian mixture model-based hidden Markov random 
field. The algorithm is implemented in MAT LAB. We also 
apply this algorithm to color image segmentation problems 
and 3D volume segmentation problems. 

1. Introduction 

Markov random fields (MRFs) have been widely used 
for computer vision problems, such as image segmenta- 
tion [ ], surface reconstruction [ ] and depth inference [ ]. 
Much of its success attributes to the efficient algorithms, 
such as Iterated Conditional Modes [^], and its consider- 
ation of both "data faithfulness" and "model smoothness" 

The HMRF-EM framework was first proposed for seg- 
mentation of brain MR images [ ] . For simplicity, we first 
assume that the image is 2D gray-level, and the intensity 
distribution of each region to be segmented follows a Gaus- 
sian distribution. Given an image Y = . . . , t/at) where 
N is the number of pixels and each yi is the gray-level in- 
tensity of a pixel, we want to infer a configuration of labels 
X = (xi , . . . , xn) where Xi e L and L is the set of all pos- 
sible labels. In a binary segmentation problem, L = {0, 1}. 
According to the MAP criterion, we seek the labeling X"^ 
which satisfies: 

X* = argmax {P(Y|X, e)P(X)}. (1) 

X 

The prior probability P(X) is a Gibbs distribution, and the 

^This work originally appears as the final project of Prof. Qiang Ji's 
course Introduction to Probabilistic Graphical Models at RPI. 



joint likelihood probability is 

p(Y|x,e) = l[p{yi\x,e) 

i 

= l[Piyi\xi,0^J, (2) 

i 

where P{yi\xi^ Ox.) is a Gaussian distribution with param- 
eters = {iixi-,crxi)' In MRF problems, people usually 
learn the parameter set 6 = {6i\l G L} from the training 
data. For example, in image segmentation problems, prior 
knowledge of the intensity distributions of the foreground 
and the background might be consistent within a dataset, es- 
pecially domain specific dataset. Thus, we can learn the pa- 
rameters from some images that are manually labeled, and 
use these parameters to run the MRF to segment the other 
images. 

The major difference between MRF and HMRF is that, 
in HMRF, the parameter set 6 is learned in an unsupervised 
manner. In a HMRF image segmentation problem, there 
is no training stage, and we assume no prior knowledge is 
known about the foreground/background intensity distribu- 
tion. Thus, a natural proposal for solving a HMRF problem 
is to use the EM algorithm, where parameter set 6 and label 
configuration X are learned alternatively. 

2. EM Algorithm for Parameters 

We still use the 2D gray-level and Gaussian distribution 
assumption. We use the EM algorithm to estimate the pa- 
rameter set 9 = {6i\l G L}. We describe the EM algorithm 
by the following [ ] : 

1. Start: Assume we have an initial parameter set 6*^^^. 

2. E-step: At the tth iteration, we have 9*^*^ and we cal- 
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culate the conditional expectation: 

Q(6|6(')) = £; [lnP(X,Y|e)|Y,e(') 



= ^P(X|Y,6(^))lnP(X,Y|6),(3) 

where x is the set of all possible configurations of la- 
bels. 

3. M-step: Now maximize (5(6 1 6*^^^ ) to obtain the next 
estimate: 

6^^+^) = argmax Q(6|6(')). (4) 



Then let 6*^^+^^ O*^*^ and repeat from the E-step. 

Let G{z;6i) denote a Gaussian distribution function with 
parameters Oi = {jii^cri)'. 



G{z;Oi 



1 / {z-jijy\ ... 
exp ( — I . (5) 
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We assume that the prior probability can be written as 



4. Calculate the posterior distribution for alU G L and all 
pixels i/i using the Bayesian rule: 



G{yi;ei)P{l\x%) 



(9) 



where is the neighborhood configuration of a;-*\ 
and 

P^'Kvi) = Y.G{yi-A)P{l\x'-Z)- (10) 

leL 

Note here we have 

P{l\xt) = |exp|- ^ Ka.xf)) .(11) 

5. Use {l\yi) to update the parameters: 



(t+i) 



(t+l)^2 _ 



i:P^'\i\yi)yi 
EP^'HiM 



(12) 



EP^'Hl\y^){yi-^^'l^'^? 
EP^^HM 



-.(13) 



P(X) = -exp(-C/(X)), 



(6) 3. MAP Estimation for Labels 



where U{x) is the prior energy function. We also assume 
that 

p(Y|x,e) = Y[p{y,\x,,e,^) 

i 



In the EM algorithm, we need to solve for X"^ that mini- 
mizes the total posterior energy 

X* = argmin {/7(Y|X, 9) + U{X)} (14) 

with given Y and 9, where the likelihood energy (also 
called unitary potential) is 



^exp(-C/(Y|X)). (7) t^(Y|X,e) = ^f/(y,|x„e) 



With these assumptions, the HMRF-EM algorithm is given 
below: 

1 . Start with initial parameter set 9*^^^ . 

2. Calculate the likelihood distribution P*^^^ {yi\xi ^Ox.). 

3. Using current parameter set 9*^^^ to estimate the labels 
by MAP estimation: 

X^^) = argmax {P(Y|X, 9^^))P(X)} 

= argmin {/7(Y|X, 9^^)) + U{X)}. (8) 

The algorithm for the MAP estimation is discussed in 
Section 3. 



E 



2(7?. 



ln( 



(15) 



The prior energy function (also called pairwise potential) 
U{X) has the form 



[/(X) = ^K(X), 



(16) 



cec 



where Vc{X) is the clique potential and C is the set of all 
possible cliques. 

In the image domain, we assume that one pixel has at 
most 4 neighbors: the pixels in its 4-neighborhood. Then 
the clique potential is defined on pairs of neighboring pix- 
els: 



(17) 
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where 



if X 2, X j 



(18) 



Note that in Eq. (17), the constant coefficient 1/2 can be 
replaced by a variable coefficient j3. We just follow [ ] to 
use the 1 /2 constant, which proves effective in many of our 
experiment results. 

We developed an iterative algorithm to solve (14): 

1. To start with, we have an initial estimate X*^^\ which 
can be from the previous loop of the EM algorithm. 

2. Provided X^^^ for all 1 < z < A/", we find 



argmin {U{yi\l) 



Vc{l,xf^)}. (19) 



3. Repeat step 2 until U{Y\X, 6) + U{X) stops changing 
significantly or a maximum k is achieved. 

4. GMM-Based HMRF 

In previous sections, we have been assuming that the in- 
tensity distribution of each region to be segmented follows 
a Gaussian distribution with parameters = {fixi^o-xj. 
However, this is a very strong hypothesis which is insuffi- 
cient to model the complexity of the intensity distribution 
of real-life objects, especially for objects with multimodal 
distributions. 

Gaussian mixture model (GMM), in contrast, is much 
more powerful for modeling the complex distributions than 
one single Gaussian distribution. A Gaussian mixture 
model with g components can be represented by parame- 
ters: 

Ol = {(/iZ,l, Cr^,l, . . . , {/ill^g, ai^g, Wl^g)}. (20) 

Compared with Eq. (5), the GMM now has a weighted 
probability 



(21) 



Now, the M-step of the EM-algorithm described in Sec- 
tion 2 changes to a Gaussian mixture model fitting problem. 
The GMM fitting problem itself can be also solved using 
an EM-algorithm. In the E-step, we determine which data 
should belong to which Gaussian component; in the M-step, 
we recompute the GMM parameters. 

5. Experiment Results 

We use the above mentioned GMM-based HMRF for 
two applications: color image segmentation and 3D vol- 
ume segmentation. For each application, minor modifica- 
tions need to be made. 





4-neighborhood in 2D 6-neighborhood in 3D 

Figure 4: Neighborhood system in 2D and 3D images. 

5.1. Color Image Segmentation 

The difference between color image segmentation and 
gray-level image segmentation is that, for a color image, the 
pixel intensity is no longer a number, but a 3 -dimensional 
vector of RGB values: Y = (y^,...,y^), and = 
{yiR^yiGjUiBy • The parameters of a Gaussian mixture 
model now becomes 

= {{^^1,1^^1,1^^1,1)^ • ' • df^l,g^^l,9^^l,9)}^ (22) 

which can be compared with Eq. (20). Also, the likelihood 
energy Eq. (15) becomes 



C/(Y|X,e) = ^f/(y,|xi,e) 
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(23) 



Example color image segmentation results are shown in 
Figure 1, 2, and 3. 

5.2. 3D Volume Segmentation 

The only difference between 2D image segmentation and 
3D image segmentation is the neighborhood system. In 2D 
images, we usually use the 4-neighborhood system or the 8- 
neighborhood system; in 3D images, we usually use the 6- 
neighborhood system or the 26-neighborhood system. The 
difference is shown in Figure 4. 

To validate our algorithm for 3D volume segmentation, 
we generate a synthetic 3D image of size 50 x 50 x 50 and 
with a foreground sphere of radius 20 at the center. The in- 
tensity of background is 0, and the foreground is 100. Ran- 
dom noise uniformly distributed within [0, 120] is added to 
the entire image, at all positions. Thus clustering methods 
such as k-means will not guarantee spatial continuousness 
of the segmentation results. A comparison of k-means seg- 
mentation and HMRF segementation is shown in Figure 5. 
With 10 EM iterations and 10 MAP iterations, while set- 
ting ^ = 1 for GMM, the 3D segmentation takes about 14 
seconds on a 2.53 GHz Intel(R) Core(TM) 15 CPU. 
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(c) 



sum of U in each EM iteration 
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ElVI iteration 



(d) 



Figure 1: Example color image segmentation results, (a) Original color image, (b) Initial segmentation by k-means. (c) Final 
segmentation by HMRF. (d) Sum of energy in each iteration. 



6. Code Documentation 

We provide the name and usage of each file of our 
MAT LAB implementation in Tabel 1 and 2, including both 
color image segmentation and 3D volume segmentation. 

7. Discussion 

In this project, we have studied the hidden Markov ran- 
dom field, and its expectation-maximization algorithm. The 
basic idea of HMRF is combining "data faithfulness" and 
"model smoothness", which is very similar to active con- 
tours [ ], gradient vector flow (GVF) [ ], graph cuts [z], 
and random walks [ ]. We also combined the HMRF-EM 
framework with Gaussian mixture models, and applied it 
to color image segmentation and 3D volume segmentation 
problems. The algorithms are implemented in MAT LAB. 
In color image segmentation experiments, we can see the 



HMRF segmentation results are much more smooth than 
the results of direct k-means clustering. In 3D volume seg- 
mentation results, the segmented object is much closer to 
the original shape than clustering. This is because Markov 
random field imposes strong spatial constraints on the seg- 
mented regions, while clustering-based segmentation only 
considers pixel/voxel intensities. 
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Figure 2: More color image segmentation results. First column: original color image; second column: initial segmentation 
by k-means; third column: final segmentation by HMRF. 
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Figure 3: More color image segmentation results (continued). 
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(a) (b) 



Figure 5: 3D volume segmentation results: (a) k-means; (b) HMRF. 



File 


Type 


Usage 


demo . m 


Runnable script 


A color image segmentation example. 
Users can run this file directly. 


image_kmeans .m 


Function 


The k-means algorithm for 2D color images. 
This will generate an initial segmentation. 


HMRF_EM.m 


Function 


The HMRF-EM algorithm. 


MRF^AP .m 


Function 


The MAP algorithm. 


gauss ianBlur . m 


Function 


Blurring an image using Gaussian kernel. 


gauss ianMask . m 


Function 


Obtaining the mask of Gaussian kernel. 


ind2i j .m 


Function 


Index to 2D image coordinates conversion. 


get_GMM.m 


Function 


Fitting Gaussian mixture model to data. 


BoundMirrorExpand . m 


Function 


Expanding an image. 


BoundMirrorShrink .m 


Function 


Shrinking an image. 


385028 . jpg 


Image 


An example input color image. 



Table 1: Name and usage of each file in color image segmentation MAT LAB code. 



File 


Type 


Usage 


demo . m 


Runnable script 


A 3D volume segmentation example. 
Users can run this file directly. 


generate_3D_image .m 


Runnable script 


Generating the synthetic input 3D image. 


image _kme an s .m 


Function 


The k-means algorithm for 3D volumes. 
This will generate an initial segmentation. 


HMRF_EM.m 


Function 


The HMRF-EM algorithm. 


MRF_MAP .m 


Function 


The MAP algorithm. 


ind2i jq . m 


Function 


Index to 3D image coordinates conversion. 


get_GMM.m 


Function 


Fitting Gaussian mixture model to data. 


Image . raw 


Raw 3D image 


An example input raw 3D image. 



Table 2: Name and usage of each file in 3D volume segmentation MAT LAB code. 
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